###################################################
### Score function for the genes on a chromosomes
###################################################

coLocScore <- function(T2,
                       locations,
                       phaseAngles,
                       bandwidth=3000000,
                       contrast=NULL){
  phaseIsTest <- FALSE
  if (class(phaseAngles)[1] == "TestResults"){
    availCont <- colnames(phaseAngles)
    if (is.null(contrast)){
      ## default to first available contrast
      contrast = availCont[1]
    }
    if (any(availCont == contrast) == FALSE){
      print("argument: contrast needs to be one of:")
      for (c in availCont){
        print(c)
      }
      stop()
    }
    else{
      phaseIsTest <- TRUE
      phaseAngels <- phaseAngles[,contrast]
    }
  }
  
  clScore <- T2
  clScore[] <- 0
  for(gene in 1:length(T2)){
    if (phaseIsTest == TRUE){
      dist <- phaseAngles[gene]
    }
    else{
      dist <- cos(phaseAngles - phaseAngles[gene])
    }
    ## browser()
    kern <- exp(-1/bandwidth^2*(locations - locations[gene])^2)
    clScoreContrib <- T2[gene]*dist*kern
    clScoreContrib[gene] <- 0 # self-score. is zero
    clScore <- clScore + clScoreContrib
  }
  return(clScore)
}


chrCoLocationTest <- function(t2, chrClassObject, phaseAngles, bandwidth){
  score <- chrClassObject@chromLocs
  useChr <- c("1","2","3","4","5","6","7","8","9","10","11","12","13","14",
              "15","16","17","18","19","X")
  for (chr in useChr){
    locations <- ChrClass@chromLocs[[chr]]
    thisChr <- intersect(names(locations), names(t2))
    if (length(thisChr) == 0){
      score[[chr]] = 0
      browser()
    }
    else{
      thisT2 <- t2[thisChr]
      thisLoc <- locations[thisChr]
      thisPhase <- phaseAngles[thisChr]
      ## browser()
      score[[chr]] <-  coLocScore(thisT2, thisLoc, thisPhase, bandwidth)
    }
  }
  return(score)
}

resampleChrCoLocatonTest <- function(t2, chrClassObject, phaseAngles,
                                     bandwidth, n){
  CoLocationScores <- chrCoLocationTest(t2, chrClassObject, phaseAngles,
                                        bandwidth)
  resampCoLocScores <- CoLocationScores
  useChr <- c("1","2","3","4","5","6","7","8","9","10","11","12","13","14",
              "15","16","17","18","19","X")
  for (chr in useChr){
    nrow <- n
    ncol <- length(resampCoLocScores[[chr]])
    resampCoLocScores[[chr]] <- matrix(0, nrow=nrow, ncol=ncol)
  }
  for(ii in 1:n){
    phaseAnglesPert <- sample(phaseAngles)
    names(phaseAnglesPert) <- sample(names(phaseAngles))
    resampCoLocScoresii <- chrCoLocationTest(t2,
                                             chrClassObject,
                                             phaseAnglesPert,
                                             bandwidth)
    for (chr in useChr){
      resampCoLocScores[[chr]][ii,] <- resampCoLocScoresii[[chr]]
    }
  }
  FN <- resampCoLocScores[["1"]]
  RN <- CoLocationScores[["1"]]
  for (chr in useChr){
    FN <- cbind(FN, resampCoLocScores[[chr]])
    RN <- c(RN, CoLocationScores[[chr]])
  }
  ## Compute Q values per gene
  qChr <- CoLocationScores
  for (chr in useChr){
    for (gene in 1:length(qChr[[chr]])){
      nn <- resampCoLocScores[[chr]] >= CoLocationScores[[chr]][gene]
      pp <- CoLocationScores[[chr]] >= CoLocationScores[[chr]][gene]
      qChr[[chr]][gene] <- median(rowSums(nn))/sum(pp)
    }
  }
  return(qChr)
}


######################################
### Plot results
######################################

plotChr <- function(qChr, ChrClass){
  require(geneplotter)
  x11()
  cPlot(ChrClass)
  useChr <- c("1","2","3","4","5","6","7","8","9","10","11","12","13","14",
              "15","16","17","18","19","X")
  for (chr in useChr){
    if(length(names(qChr[[chr]])[qChr[[chr]] < 0.01]) > 1){
      cColor(names(qChr[[chr]])[qChr[[chr]] < 0.01], "red", ChrClass)
    }
  }
}
